
clear
do "...\First.do"

********************************************************************************

clear
use "$work\full_sample.dta"

keep if inrange(age,40,70)
keep if inrange(year,1995,2019)
sort pnr year

********************************************************************************
destring udd_level, g(udd_kode)
tab udd_kode

g udd_group=1 if inlist(udd_kode,10,15) 
replace udd_group=2 if inlist(udd_kode,20,25,35)
replace udd_group=3 if inlist(udd_kode,40,50)
replace udd_group=4 if udd_kode>50
replace udd_group=. if udd_kode==.

lab define udd_group2 1 "Primary education" 2 "Vocational Education" 3 "College" 4 "University"
lab values udd_group udd_group2


rename death_cancer_lung death_lung

rename GP_visit GP_contact
rename N_visits2 N_visits
rename services mean_services
rename COPD_med meds_COPD

foreach var of varlist death death_CVC death_cancer death_lung  married yder_ses_max GP_contact N_visits mean_services total_cost first_lung  statins ACE diabetes_control metformin meds_COPD ACSC_COPD  {
	
reg `var' i.age male i.year
predict res_`var', res

bys low_ses: egen mean_res_`var'=mean(res_`var')
egen base_`var'=mean(`var')

g controlled_`var'=mean_res_`var'+base_`var'

replace controlled_`var'=controlled_`var'*100
bys low_ses: egen mean_`var'_controlled=mean( controlled_`var')

bys low_ses: egen mean_`var'=mean(`var')

drop mean_res_`var' base_`var'

}


keep pnr male yob dk low_ses mean_* mean_*_controlled mean_*
drop mean_age mean_dk mean_grades mean_services mean_male

duplicates drop

save "$work\temp_sumstat_fig2.dta", replace


********************************************************************************
clear
use "$work\temp_sumstat_fig2.dta"


foreach var of newlist  death death_CVC death_cancer death_lung  married yder_ses_max GP_contact N_visits mean_services total_cost_GP first_lung statins ACE diabetes_control metformin meds_COPD ACSC_COPD  {
ttest  mean_`var'_controlled, by(low_ses)
g mean0_`var'_c= round((`r(mu_2)'-`r(mu_1)')*1000/1000, 0.00001)
g base0_`var'_c= round((`r(mu_1)')*1000/1000, 0.00001)
g base1_`var'_c= round((`r(mu_2)')*1000/1000, 0.00001)

}


keep mean_* mean0_* base0_* base1_* low_ses
drop mean_* low_ses
duplicates drop

rename *_c *

g id=1

reshape long mean0_ base0_ base1_, i(id) j(var) string

drop if var=="yder_ses_max"
set obs 24
replace var = "Mortality" in 24

set obs 25
replace var = "Care" in 25

set obs 26
replace var = "Reimbursement" in 26

set obs 27
replace var = "Health behavior" in 27


replace var="All Cause Mortality" if var=="death"
replace var="COPD" if var=="death_COPD"
replace var="CVC" if var=="death_CVC"
replace var="Diabetes" if var=="death_diabetes"
replace var="Cancer" if var=="death_cancer"
replace var="Lung cancer" if var=="death_lung"

cap drop sort
g sort=1 if var=="Mortality"
replace sort=2 if var=="All Cause Mortality"
replace sort=3 if var=="CVC"
replace sort=4 if var=="Cancer"
replace sort=5 if var=="Lung cancer"


replace sort=8 if var=="Care"
replace sort=9 if var=="GP_contact"
replace sort=10 if var=="N_visits"
replace sort=11 if var=="mean_services"
replace sort=14 if var=="total_cost_GP"

replace var="PCP visit (D)" if var=="GP_contact"
replace var="PCP visits (N)" if var=="N_visits"
replace var="Services per visit" if var=="mean_services"


replace var="Total reimbursement" if var=="total_cost"
replace var="PCP reimbursement" if var=="total_cost_GP"
replace var="Specialist reimbursement" if var=="cost_special" 


replace sort=17 if var=="Health behavior"

replace sort=20 if var=="statins"
replace sort=21 if var=="ACE"
replace sort=23 if var=="meds_COPD"
replace sort=24 if var=="ACSC_COPD"
replace sort=25 if var=="metformin"
replace sort=26 if var=="diabetes_control"
replace sort=27 if var=="first_lung"


replace var="Statins" if var=="statins"
replace var="COPD medication" if var=="meds_COPD"
replace var="ACSC COPD" if var=="ACSC_COPD"
replace var="Metformin" if var=="metformin"
replace var="Diabetes check-up" if var=="diabetes_control"
replace var="Lung scan" if var=="first_lung"

cap drop procent
g procent=mean0_/base0_*100
g sort2=_n
sort sort

insobs 1, before(6)
insobs 1, before(12)
insobs 1, before(1)

replace sort=1 if sort[_n+1]==1
replace sort=8 if sort[_n+1]==8
replace sort=16 if sort[_n+1]==17


graph hbar (mean) procent if sort!=., over(var, sort(sort) label(labsize(small))) graphregion(color(white) ilcolor(white) lcolor(white)) bgcolor(white) ytitle("SES gradient %", size(small)) bargap(1.5) blabel(bar, format(%6.1f)) ylabel(-0(20)100, labsize(small)) missing bar(1, color(sienna%40)) ///
legend(order(1 "{it: {Mortality}}" 2 "All cause" 3 "CVC"))


graph export "$fig\gradient_graph_new_adjusted.png", replace 
graph export "$fig\gradient_graph_new_adjusted.pdf", replace 

